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The linking number (topological entanglement) and the writhe (geometrical entanglement) of a 
model of circular double stranded DNA undergoing a thermal denaturation transition are investi- 
gated by Monte Carlo simulations. By allowing the linking number to fluctuate freely in equilibrium 
we see that the linking probability undergoes an abrupt variation (first-order) at the denaturation 
transition, and stays close to 1 in the whole native phase. The average linking number is almost zero 
in the denatured phase and grows as the square root of the chain length, N, in the native phase. 
The writhe of the two strands grows as \/7V in both phases. 



I. INTRODUCTION 

It is well known that topological (knotting and link- 
ing) and geometrical (supercoiling) entanglements in long 
DNA molecules are critical to the functioning of the 
cell 1^. For this reason there exist enzymes such as 
topoisomerases and recombinases which facilitate cellular 
metabolism by changing the geometry and the topology 
of DNA |-|. 

For closed loops of double-stranded DNA (dsDNA) the 
consequences of the double-helix linking number Lk being 
non zero have been extensively studied |^-^. For exam- 
ple, it is known that a large imposed linking number leads 
the chain to twist upon itself. This is the phenomenon of 
supercoiling in DNA, analogous to the familiar buckling 
of a twisted tube or wire . The linking entanglement 
is believed to play an important role in structural tran- 
sitions of double stranded chains such as local denatu- 
ration and cruciform structures formation. In particular 
at the melting transition it is possible that progressive 
supercoiling of the native part of the molecule affects the 
denaturation process and there are indeed experimental 
indications that supercoiled structures denature at higher 
temperature and over a broader temperature range than 
DNA molecules in the relaxed state ||ll|,|l2|. However, 
despite this evidence of the effects that topology and ge- 
ometry can have on the denaturation transitions of circu- 
lar dsDNA's, theoretical studies of thermal denaturation 
have in general neglected this aspect up to know. Indeed, 
starting with the seminal work of Poland and Sheraga 
[^3|, statistical mechanics studies of denaturation have 
focussed mainly on the nature of the transition and on 
how this depends on properties such as the self and mu- 
tual avoidance of the strands 1 14 l6| , the bending rigidity 
of the bound segments |0] , or the inhomogeneity of the 



base pair sequence An attempt to include topologi- 
cal constraints in statistical models of thermal denatura- 
tion has been made recently by Rudnick and Bruinsma 
|l9| . This study, carried on within the Poland and Sher- 
aga model, is based on the introduction at mean field 
level of an elastic strain energy concentrated in the native 
part of the chain. Their model completely neglects the 
3D nature of the entanglement complexity of the double 
stranded chain. In order to reach more firm conclusions 
about the effects of the torsional strain on the nature 
of thermal denaturation, one needs of course models in 
which these effects coexist with a realistic representation 
of the geometrical and topological entanglement of DNA. 

In the present article we address the problem of de- 
scribing properly the topological and geometrical com- 
plexity of DNA when it undergoes the denaturation tran- 
sition. More precisely, we set up a model in which, while 
denaturation is induced by energetic factors not directly 
related to the entanglement complexity of the chain, the 
topological and geometrical properties can be meaning- 
fully represented and monitored in equilibrium. This ap- 
proach, which has some complementarity with respect to 
that of Ref. , should be seen as a first step towards de- 
scriptions that embody also torsional strain interactions 
induced by a linking constraint, which are not considered 
here. 



II. MODEL AND PHASE DIAGRAM 

We consider a model of circular double stranded 
DNA where the two strands are mimicked by two non- 
overlapping self-avoiding polygons (SAP) of N edges on 
the cubic lattice. Suitable short range attractive interac- 
tions between homologous pairs of monomers of the two 
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strands are then considered to induce the formation of 
the double structure. Each SAP is rooted at one ver- 
tex and the rooted vertices of the two SAP's are fixed 
to remain always a lattice distance apart. We identify 
the configuration uj of the two SAP's by the positions 
ri{i) and ^2(1), i = 1, 2, . . . A^, of the visited lattice ver- 
tices. A binding energy — e, (e > 0) is associated to 
each pair of vertices (ri(i), r2(j)) having i = j and dis- 
tance \ri{i) — r2{j)\ < \/3 in lattice units. We call such 
kind of pair a contact. The extension of the range of 
the interaction to the third neighbors is suggested by the 
need of conferring a reasonable degree of flexibility to 
the two SAP's when they are bound together to form 
the double chain structure. At low temperature, shorter 
ranges of interaction would indeed cause double stranded 
(ribbon-like) structures that are quite rigid and, even 
if rigidity effects are not expected to be relevant for the 
asymptotic features of the denaturation transition [ p^ , 
their reduction can avoid slow crossover effects. 
The Hamiltonian of a dsDNA configuration uj is 



H{ll!) = —771(10) 



(1) 



where 171(0;) is the number of contacts as defined above. 
The canonical partition function Z at an inverse temper- 
ature = 1/r is 



(2) 



and most thermodynamic quantities follow as weighted 
averages normalized by Z . The binding energy is taken to 
be the same all along the strand, i.e., we neglect the het- 
erogeneity of base pair interactions of specific sequences. 
This is a reasonable first approximation if we consider 
that a single unit lattice step should correspond to a per- 
sistence length of the strands. So e represents in fact 
an average binding energy for a set of several base pairs. 
At = the two strands behave as independent SAP's, 
whereas at sufficiently high values of (3 the two SAP's are 
bound, i.e a macroscopic number (cx N) of contacts oc- 
curs. In analogy with similar previously studied models 
| p^ , pO[ , a melting transition is expected to occur at some 
f3 = Pc- Coming from low temperatures, this transition is 
driven by the formation of denatured loops, whose length 
can be measured by the number of unbound monomers t 
comprised between two consecutive contacts of the cor- 
responding strand segments ||l^,|l6| . 

Configurations have been sampled by Monte Carlo 
(MC) methods. Since the system is strongly interact- 
ing we have adopted a multiple Markov chain approach 
in which one samples simultaneously at a variety of dif- 
ferent temperatures and "swaps" configurations between 
contiguous temperatures. The swap probability is cho- 
sen so that the limit distribution of the process is the 
product of the Boltzmann distributions at the individual 
temperatures [ pT|] . The underlying (symmetric) Markov 
chain used is based on a combination of pivot moves for 
SAP p2] and of local moves [ESl . In addition, to increase 



the mobility of the MC sampling in the native phase, 
we have introduced a new, "double inversion" move (see 
Fig. |l|). This move has the advantage that it may change 
the mutual entanglement of the two SAP's while keeping 
the number of contacts fixed; it is thus particularly effec- 
tive in sampling configurations in which the two strands 
are strongly bound and linked. The result of applying 
this move is a net increase of the mobility of the MC in 
the low temperature, native phase. 



FIG. 1. Schematic representation of the double-inversion 
move; two monomers i and j are selected randomly (with 
i < j) along the SAP sequence and a move is tried if the 
positions f\(i) and r\(j) (first SAP) and r2(i) and r2(j) (sec- 
ond SAP) satisfy the relation f\(i) — r2(i) = ri(j) — r2(j). 
The proposed move is an inversion of the portion of the two 
chains between i and j with respect to their middle point 
f™ = l/2(ri(i)+f2(j)). 

In our simulations, each MC step consists of 0(1) pivot 
and double inversion moves and 0(N) local moves. By 
running in parallel a number M « 30 of Markov chains 
we are able to obtain good sampling for chains up to 
2N = 1200, from /3 = up to /3 = 1.1. For all lengths we 
obtain samples with at least 10^ independent data points 
for almost all f3's considered. Only for the highest /3's the 
correlation of successive samples reduces the number of 
independent data points to ^ 10^. To interpolate the ob- 
tained data at intermediate temperatures we have used 
the multiple histogram method |Q . 

In order to study the entanglement of the configura- 
tions as a function of temperature and in particular at 
the melting transition, we need first to characterize the 
phase diagram of the model and to estimate Pc- Two 
regimes are manifestly present: a low /3 regime charac- 
terized by two unbound and independent SAP's, and a 
high P phase dominated by configurations in which the 
two strands are strongly bound (see Fig. ^). 

A transition between these two regimes is suggested by 
the behavior of the average energy per monomer 
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N ZN ^ ^ ' 



(3) 



as a function of (3 (see Fig. 
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FIG. 4. Plot of the specific heat C as a function of /3, for the 
same values of N as in Fig. ^. Note that each curve displays 
a maximum that increases as A'^ increases. Since we expect 
the transition to be asymptotically first-order (0 = 1, c > 2), 
the peaks should develop into a 5- like singularity as N oo. 




Indeed from Fig. || we note that, for P < Pc ~ 0.9, 
{H)/N goes to zero as ^ cxd whereas {H)/N reaches 
a finite limit for (3 > (3c- The crossover between these 
regimes is particularly sharp and this is confirmed by the 
behavior of the specific heat 



^ 13^ d{H) 
~ N dl3 



(4) 



as a function of /3, which is plotted in Fig. |4|. In the 
proximity of /3c the singular part of the specific heat is 
expected to scale as 



C^N^<^-^!{{(3-(3c)N*) 



(5) 



FIG. 2. Snapshots of typical configurations at /3 lower (top) 
and higher (bottom) than the denaturation inverse tempera- 
ture fic- 




FIG. 3. Plot of the energy density {H)/N as a function of 
13, for N = fOO(o), 150(x), 200(a), 300(D), 400(V), 500(o), 
600(»). As A'^ increases, two behaviors are clearly distin- 
guished: for /3 < 0.9, {H)/N 0, while {H)/N ^const 
for /3 > 0.9. 



where / is a scaling function and 4> is the crossover 
exponent. Since for large N one has Craax{N) = 
max/3 C(Af, /3) ^ N^'l'^^, by a linear fit of InCmax vs 
InA^ we obtain (p — 0.9(1). From Eq. (||) we expect 
also {f3 — Pc) ^ N~'^ and by extrapolating Pc{N) vs 
l/N"^, with (j) given by the previous estimate, we obtain 
/3c = 0.905(5). 

As in other recently studied models of DNA denatu- 
ration 1 15 1^, it is useful to look at the probability dis- 
tribution P{£, N) of the denatured loops of length £. In- 
deed, it turns out that for /3 < /3c Pie,N) - e^'=f{l/N) 
with an exponent c connected to (j) by the relation (p = 
min{c — 1,1}, if c > 1. Hence, a continuous transition 
corresponds to c < 2, while c > 2 gives ip — 1 and a 
first order denaturation. To estimate c, we examine the 
P{£,N = 600), shown in Fig. |, for different /3's. Its 
slope at /3 = 0.9 and in the range 10 < £ < 180 has 
the borderline value c = 2.01(5). This value is slightly 
lower than the one predicted for somewhat simpler mod- 
els considered recently, in which the geometry of contacts 
is drastically simplified to the extent that linking entan- 
glement cannot be defined for the two strands |lq . This 



3 



small discrepancy is most probably due to the relatively 
shorter chains considered here and to the less accurate 
determination of /3c- We believe that the transition has 
a first order character. This conclusion is also supported 
by the behavior of the linking probability discussed be- 
low. It should also be taken into account that the present 
model has interactions which induce some bending rigid- 
ity of the double stranded structure. These effects have 
already been shown in to produce a slight lowering 
of the effective c for finite N. 
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FIG. 5. A log-log plot of the loop probability distribution 
P{£,N = 600), for different f3 values. For /3 < 0.9 the curves 
display a power law behavior, while for (3 = 1.06 strong cut-off 
effects are present. 



III. MUTUAL ENTANGLEMENT: LINKING 
NUMBER 

The topological and geometrical properties of the 
model are then studied by looking at the behavior, as 
a function of T, respectively of the linking probability 
and of the writhe, which gives a measure of the degree of 
supercoihng pst . 

In this section we ask for the probability that the two 
strands are linked as a function of T and, in particular, 
its behavior at the melting transition, whose location has 
been estimated in the previous section. First let us be 
clear on what we mean by two linked curves. 

Two disjoint simple closed curves Ci andC2 are topolog- 
ically unlinked if there is a homeomorphism of onto 
itself, H : B? ^ B? , such that the images H{Ci) and 
H{C2) can be separated by a plane |^. 

A weaker notion of linking is given by the homological 
linking definition. Ci is homologically unlinked from C2 if 
Ci bounds an orientable surface which is disjoint from C2- 
Homological linking is a symmetric relation, and homo- 
logical linking implies topological linking. In this paper 
we shall be concerned just with homological linking. The 
reasons to restrict ourselves to this definition of linking 
are twofold. In first place this is the simplest topologi- 
cal property to be checked computationally. In addition, 
despite the fact that homological linking is the weakest 
form of linking, it is known to be a good indicator of 



topological linking for configurations in which the two 
SAP's are strongly interpenetrating |27|. This is the sit- 
uation that we expect to occur in the proximity of the 
denaturation transition. 

A method to detect whether or not the two strands are 
homologically linked consists in orienting each of the two 
strands Ci and C2, and to project them onto a plane so 
that no vertex in the projection of Ci coincides with any 
vertex in the projection of C2 , or viceversa (regular pro- 
jection). In this way the mutual crossings are transverse 
(no edge overlaps between the two projections) and at 
each point where Ci crosses under C2 we assign a value 
-1-1 or —1, according to the orientation of the crossing 
(see Fig |). The sum of these crossing numbers is called 
the linking number of the two curves , Lk(Ci,C2), and 
the two curves are homologically linked if and only if 
Lk{Ci,C2)^0^. 

FIG. 6. Positive and negative crossings determined a by 
left-hand rule. 




0.2 0.4 0.6 0.8 



FIG. 7. Plot of the linking probability as a function of 0. 
Different curves correspond to different A'^ values, as in Fig. H. 
In correspondence to the transition {/3 ~ 0.9) an abrupt jump 
of Phom develops. 

In Fig. 1^ we show the probability, Phom, for the two 
circular strands to be homologically linked (i.e. having 
non-zero linking number) as a function of (3 and for dif- 
ferent N values. The qualitative trends are clear. The 
linking probability increases with increasing (3 at fixed 
N. In particular at high /3, (i.e. in the native phase), the 
linking probability is very close to 1, indicating that the 
double stranded structures are very likely to be linked 
even for finite values of N. It is interesting to notice 
that, in the proximity of the denaturation transition, the 
change in slope is very sharp, suggesting a possible dis- 
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continuity in the value of Phom- This possible feature is 
more evident in Fig. || where Phom has been plotted as 
a function of for different /3 values. It seems clear 
that for (3 < (3c i.e. above the denaturation transition, 
the curves extrapolate (as N oo) to values p = p{P) 
with p varying gradually, but sensibly, with /?, whereas as 
soon as (3 > /3c (i.e. in the bounded region) the limiting 
linking probability is very close to 1. At /3 = the linking 
probability seems to extrapolate to a value p{0) ~ 0.2. 
A more detailed analysis of the linking probability can 
be carried out by looking at the probability, PhkiP), to 
have a given linking number Lk. In Fig. PLki/3) has 
been plotted as a function of /? for different Lk values 
and for N — 400. Again one can observe that at /3 ~ /3c 
there is an abrupt decay of Po{f3) (i.e. the probability 
to be unlinked) to which corresponds a sharp increase 
in all the probabilities of having nonzero values of Lk. 
For /3 >> /3c it seems that Pa{(3) goes to a constant dif- 
ferent from zero; this is however a finite size effect (the 
plot is for N — 600) and as N increases Pq{P) decreases 
monotonically. 
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FIG. 8. Plot of Phom as a function of 1/N for different (3 
values. By observing that (3 = 0.88 is below the transition 
point /3c = 0.905(5), while (3 — 0.91 is just above /3c, it is ev- 
ident that a change in the asymptotic behavior of the curves 
takes place in proximity of the transition. 
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FIG. 9. Plot of the probability of having a fixed link- 
ing number Lfc > as a function of /3 for TV = 400. In 
the plot we have just considered positive values of Lk since 

P^Lk{l3) = PLk{l3). 
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FIG. 10. Average of the absolute value of Lk as a function 
of (3. Difi'erent curves correspond to different A*' values, as in 
Fig. ^ Notice that for (3 > 13c there is an abrupt increase of 
the average linking number. The scaling of is consis- 

tent with the form {\Lk\) ~ VN and with /3c = 0.905(5), as 
shown in the inset. 

In Fig. |l^ we show the average of the absolute value of 
the linking number (|iA;|) as a function of (3 for various 
values of A^. In the whole range of /3 < /3c and for each 
A^, (|i^|) is very close to zero and grows very slowly with 
A^. As /3 > /3c there is an abrupt change to a regime char- 
acterized by a more rapid increase ol {\Lk\) as a function 
of A^. By assuming that for (3 > Pc the following power 
law behavior holds 



\Lk\ 



(6) 



The evidence of a discontinuity of Phom presented 
above supports the first order character of the denatu- 
ration transition of the model. This should be expected 
in view of the results of Refs. ||l5|,|l6| and of the fact that 
Lk is unrestricted and no contributions to the energy are 
associated here to the twist of the bounded DNA seg- 
ments. 



we can estimate the a exponent through a simple linear 
fit of log(|Lfc|) vs logA^. This gives 



cr = 0.52 ± 0.01 



(7) 



for the whole range /3 > /3c (see Fig. 11). It is inter- 
esting to notice that such estimate is in good agreement 
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with the one obtained for the two boundary curves of an 
orientable lattice ribbon model ||2^. In other words the 
low temperature phase of our model presents the same 
asymptotic topological properties as those of a lattice 
ribbon model introduced some time ago to describe the 
entanglement complexity of double-stranded molecules in 
a good solvent p9| , ^ . This is what one should expect 
for a model of dsDNA in the native state and confirms 
the adequacy of our model. 

A further indication of the scaling behavior in Eq. (|^) 
is shown in the inset of Fig. where we plot {\Lk\)/\fN 
vs /3: clearly in the range (3 > (3c all the data collapse 
onto a single curve whereas in the denatured phase (i.e. 
(3 < Pc) the curve approaches zero as N increases. 



IV. GEOMETRICAL ENTANGLEMENT: 
WRITHE 

In this section we analyze the geometrical entangle- 
ment of the system as a function of temperature by com- 
puting the writhe of each SAP separately. In order to de- 
fine a writhe, consider any simple closed curve in R? ^ and 
project it onto B? in some chosen direction. In general 
the projection will have crossings and, for almost all pro- 
jection directions, these crossings will be transverse, so 
that we can associate a sign or —1 with each crossing. 
For this projection we perform the sum of these signed 
crossing numbers and average over all possible projection 
directions. This average quantity is the writhe Wr of the 
curve 0. The writhe is a geometrical quantity (since it 
is not invariant under ambient isotopy) and, contrary to 
the integer Lk^ it is a real number measuring the extent 
to which the strand is supercoiled. In principle one needs 
to average the sum of the signed crossing numbers over 
all (regular) projections but there is a useful theorem ||3l| ] 
which considerably simplifies the calculation for polygons 
on the simple cubic lattice. Indeed for a SAP in the cubic 
lattice, the computation of the writhe equals the average 
of the linking numbers of the given SAP and its push-offs 
(translate through a sufficiently small distance) into four 
non-antipodal octants 

For a single non-interacting SAP it is known that the 
mean of the absolute value of the writhe behaves as 



{\Wr\) - N^, 



(8) 



and we expect the same power law behavior for our model 
in the denatured regime (i.e. j3 < (3c) where the two 
strands behave essentially as two independent SAP's. As 
a check we estimate the C exponent at (3 ~ Q assuming 
the scaling behavior (H). By a linear fit of the log- log 
plot of the N dependence of (| Wr\) at /3 = 0, and with 
N e {300, 400, 500, 600}, we obtain 

C = 0.505(4), (9) 

in full agreement with the corresponding estimate for a 
single SAP, which is also known to be bounded from be- 
low by 1/2 [Hj. Next, by assuming the validity of (||) 



at any /3, we estimate the exponent C as a function of 
(3 (see Fig. |ll|). With the exception of a little variation 
in the proximity of the transition, it seems clear that C 
remains very close to the (3 — estimate for all (3 values. 
Thus, the denaturation transition does not seem to af- 
fect the exponent of the power law dependence on N of 
the absolute value of the writhe. This is not surprising 
if one observes that for (3 < (3c the two circular strands 
behave essentially as two independent SAP's, whereas 
for (3 > (3c the two strands are tightly bound together 
forming a ribbon-like structure. Indeed the square root 
of N behavior of the absolute value of the writhe has 
been found also for the two boundaries of an orientable 
lattice ribbon model p8|-^. This confirms, also from 
the geometrical point of view, the good correspondence 
between the native regime of our model and the lattice 
ribbon model for double stranded molecules. The devia- 
tions of C from 1/2 observed in the neighborhood of (3c are 
probably due to finite size effects. However, we can not 
exclude that a peculiar, distinct scaling regime prevails 
right at the transition. We know that there the statistics 
of denaturated loops is completely different from that at 
(3 > (3c and this could imply a modification in the scaling 
behavior of (| Wr\). To establish whether a distinct scal- 
ing of (I Wr\) exists would be, however, very challenging 
numerically and beyond the scope of the present work. 
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FIG. 11. Effective exponent a of the linking number (A) 
and C, of the writhe (l>) as a function of (3. The variations 
of close to the transition are probably due to corrections to 
the scaling Eq. that are much larger than the statistical 
errors displayed in the figure. 



V. CONCLUSIONS 



In this paper we have studied the topological and geo- 
metrical entanglement of a lattice model for circular dou- 
ble stranded DNA undergoing a denaturation transition. 
We have shown that, in the limit of very long chains, the 
linking probability between the two strands is a function 
of (3 in the denaturated phase, whereas it jumps very 
rapidly to values close to unity as soon as /3 > /3c, i.e. 
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in the bound state. This feature is confirmed by the be- 
havior of the average hnking number in the two phases: 
it turns out that is a small constant for (3 < (3c 

and grows roughly as \/7V in the low T phase. This sug- 
gests that our model in the native phase is, as far as 
the homological linking is concerned, similar to a ribbon 
model. This analogy is also confirmed by the behavior 
of the absolute value of the writhe as a function of (3. 
In particular (| Wr|) roughly scales with N as N^-^^ for 
all /3's and no supercoiling effects have been observed in 
this model of denaturation. This can be explained as fol- 
lows: in our model the linking number is not fixed to a 
particular value and can fluctuate freely in the equilib- 
rium statistics. Experimentally this would correspond, 
at least qualitatively, to the presence in the solution of 
the topoisomerases whose function is to change the link- 
ing number continuously. A different scenario can show 
up if instead the system is constrained to have a fixed 
linking number (no topoisomerases in solution). Indeed 
the fact that, for our model with unconstrained topol- 
ogy, the probability of being linked undergoes, right at 
denaturation, an abrupt (first order like) jump, suggests 
that the imposition of a constraint on the topology of 
the dsDNA molecule (by fixing the linking number be- 
tween the two SAP's in our model) would affect rather 
sensibly some features of the transition, such as the melt- 
ing temperature. There are indeed experimental indica- 
tions that supercoiled structures characterized by a large 
(fixed) linking number, display a melting transition at 
higher temperatures than DNA molecules in the relaxed 
state fill. 
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